% fig2 相图(Phase)

clear; clc; close all;

%% —— 参数网格 & 初始设置 —— 
g12 = 1;   g21 = -2;  g22 = 2;
mu_vec  = linspace(0,1,500);    % μ 网格
g11_vec = linspace(-4,4,500);   % g11 网格
nMu  = numel(mu_vec);
nG11 = numel(g11_vec);

Niter  = 600;   % 总迭代步数
Ntrans = 400;   % 热身步数（丢弃）

x0 = 0.1;  y0 = -0.1;            % 轨道初值
tolLE   = 1e-4;                  % 判断 LE≈0 的阈值
tolUniq = 1e-4;                  % 周期点唯一化精度

% 预分配：存放每个 (g11,μ) 的状态编码
state_map = zeros(nG11, nMu);

%% —— 并行扫参，内联计算前两大 LE 和周期数 —— 
parfor ii = 1:nG11
    g11 = g11_vec(ii);
    row = zeros(1, nMu);
    for jj = 1:nMu
        mu = mu_vec(jj);
        % 初始化轨道与扰动向量
        x = x0;  y = y0;
        v1 = [1;0];  v2 = [0;1];
        s1 = 0;     s2 = 0;
        y_traj = zeros(Niter - Ntrans,1);
        idx_y = 0;
        % 迭代与 LE 计算
        for k = 1:Niter
            x_next = mu*x + g11*sin(x) + g12*sin(y);
            y_next = mu*y + g21*sin(x) + g22*sin(y);
            J = [mu+g11*cos(x), g12*cos(y);
                 g21*cos(x),    mu+g22*cos(y)];
            % Gram–Schmidt 正交化
            w1 = J*v1; n1 = norm(w1); v1 = w1/n1;
            w2 = J*v2 - (v1'*(J*v2))*v1; n2 = norm(w2); v2 = w2/n2;
            if k > Ntrans
                idx_y = idx_y + 1;
                y_traj(idx_y) = y_next;
                s1 = s1 + log(n1);
                s2 = s2 + log(n2);
            end
            x = x_next;  y = y_next;
        end
        % 计算 LE1, LE2
        M = Niter - Ntrans;
        LE1 = s1 / M;
        LE2 = s2 / M;
        % 周期检测：唯一化 y_traj
        pts = unique(round(y_traj/tolUniq)*tolUniq);
        pnum = numel(pts);
        % 状态分类
        if LE2 > 0
            code = 2;            % HC
        elseif LE1 > tolLE
            code = 1;            % CH
        elseif abs(LE1) < tolLE && pnum > 8
            code = 0;            % QP
        else
            if pnum == 1
                code = 11;       % SP
            elseif pnum>=2 && pnum<=8
                code = 10 + pnum; % P2→12, …, P8→18
            else
                code = 20;       % MP
            end
        end
        row(jj) = code;
    end
    state_map(ii,:) = row;
end

%% —— 重映射编码到 1…12 索引 —— 
codes = [0,1,2,11,12,13,14,15,16,17,18,20];  % 原始编码
idx_map = zeros(size(state_map));
for k = 1:numel(codes)
    idx_map(state_map==codes(k)) = k;
end

%% —— colormap —— 
cmap = [
    1.00, 0.50, 0.00;  %  1: QP  橙
    1.00, 0.00, 0.00;  %  2: CH  红
    1.00, 0.75, 0.80;  %  3: HC  浅粉
    0.00, 0.00, 0.00;  %  4: SP  黑
    0.00, 0.00, 0.50;  %  5: P2  深蓝
    0.00, 0.00, 1.00;  %  6: P3  蓝
    0.25, 1.00, 0.25;  %  7: P4  黄绿
    1.00, 1.00, 0.00;  %  8: P5  黄
    1.00, 0.50, 1.00;  %  9: P6  品红
    0.00, 1.00, 0.00;  % 10: P7  绿
    0.50, 0.50, 0.50;  % 11: P8  灰
    0.00, 1.00, 1.00;  % 12: MP  青
];

%% —— 绘图 —— 
figure;
imagesc(mu_vec, g11_vec, idx_map);
set(gca,'YDir','normal');

colormap(cmap);
caxis([1 12]);
hcb = colorbar(...
    'Ticks',1:12, ...
    'TickLabels',{'QP','CH','HC','SP','P2','P3','P4','P5','P6','P7','P8','MP'}, ...
    'Location','eastoutside');
hcb.FontSize = 8;

xlabel('\mu','FontSize',14);
ylabel('g_{11}','FontSize',14);
title('Fig.2 完整状态分类 (g_{12}=1)','FontSize',16);

labels = {'QP','CH','HC','SP','P2','P3','P4','P5','P6','P7','P8','MP'};

total = numel(state_map);
for k = 1:numel(codes)
    cnt = sum(state_map(:) == codes(k));
    pct = cnt/total*100;
    fprintf('%s: %d points, %.2f%%\n', labels{k}, cnt, pct);
end